C- f for sm calls on HPUX and AIX, use '' for others
      Program rho0Yudinnmain
C- _TRACE "@wterm' IB=',IB,' Ncheb=',Ncheb,' nDersMAX=',nDersMAX," ;
C- uses eosyudin.trf
C- constructs stationary barotropic configurations in fast rotation along rho0 series
C- this code may be used for testing non-rotating models
C- in this test one can put delta==1, while normally for alpha-law delta==2
C- but now delta=2 is not plotted correctly
C- the method is described in AA290(1994)674 for ursos
C- main variables:
C- phi(*,*) - grav.potential \hat\Phi as in eq.2 and 10
C- Fpsi - centrifugal potential \hat\Psi as in eq.8 and 14
C-*****************************************************************
      use utilities_module
C      use,intrinsic :: iso_c_binding
      use dflib
      use cylinderGrid
      use accur
C-        _include accur;
C-        include 'accur.for';
C-        double precision  pi;
C:  for some functions *
      interface
      doubleprecision Function Fpsi(rcylw,rotlaw)
      implicit none
      doubleprecision::rcylw
      character*(*)::rotlaw
      endfunction
      doubleprecision Function Fdpsi_drcyl(rcylw,rotlaw)
      implicit none
      doubleprecision::rcylw
      character*(*)::rotlaw
      endfunction
      doubleprecision Function chi(rcylw,rotlaw)
C- \hat\Psi^0
      implicit none
      doubleprecision::rcylw
      character*(*)::rotlaw
      endfunction
      doubleprecision Function powrsmart(x,alpha)
      implicit none
      doubleprecision::alpha,x
      endfunction
      endinterface
C-       include 'statd.for';
C-        _include grid;
      doubleprecision S_mass,G_N,urho,uH
      Parameter(S_mass=1.989d33,G_N=6.6732d-8)
      Parameter(urho=2d9,uH=1.7d18)
C- WD
C: parameters like S_mass, Imax0 and commons *
C-       Parameter(S_mass=1.989d33, G_N=6.6732d-8); -- now in WDunits.inc
C       integer NrDAT,NzDAT,NtDAT,indmin(1),indmax(1);
C       Parameter(NrDAT=IA,NzDAT=4*IA,NtDAT=JMax);
C       REAL X(NrDAT+1),Y(NzDAT+1),phiRcylZplot(NrDAT,NzDAT);
C--    NrDAT=IA;   -- points on Rcylindrical
C--    NzDAT=4*IA; -- points on Z -- more points better for flattened bodies
C--    NtDAT=JMax; -- points on thet
C       integer NDIMP,NDP,NIP,ir,jt,i,j, -- for spherical coords
C               icyl,jcyl; -- for cylindr. coords
C       Parameter(irho0max=@irho0max);
C       double precision  rspl(0:NrDAT),thet(0:NtDAT),Phirthet(0:NrDAT,0:NtDAT),
C       rhoRcylZ(0:NzDAT,0:NrDAT)/,rhodV(0:NzDAT)*
C-     E01DAF  based on Example Program Text
C-     Mark 14 Release.  NAG Copyright 1989.
C-     .. Parameters ..
C-      INTEGER          NIN, *;
C-       PARAMETER        (NIN=5,*=6);
      INTEGER M,LIWRK,LWRK
      PARAMETER(LIWRK=NrDAT+1+2*(NrDAT-2)*(NtDAT-2),
     *	LWRK=(NrDAT+7)*(NtDAT+7))
      INTEGER ChMMAX,ChNMAX,nDerMX,ChLWRK,ChLIWRK
      PARAMETER(ChMMAX=NrDAT+1,ChNMAX=NrDAT+1,
C- ChNMAX==ChMMAX when no derivatives are given
     *nDerMX=2,ChLWRK=7*ChNMAX+5*nDerMX+ChMMAX+7,ChLIWRK=2*ChMMAX+2)
      INTEGER Ncheb,nDers(ChMMAX),skipCheb
      DOUBLEPRECISION Xch(ChMMAX),Ych(ChNMAX)
      DOUBLEPRECISION chebA(ChNMAX),ChWRK(ChLWRK),resCheb(ChMMAX),
     *	pt(ChMMAX)
C-     .. Local Scalars ..
C-       DOUBLE PRECISION STEP, XHI, XLO, YHI, YLO;
      INTEGER IFAIL,MX,MY,PX,PY
C-     .. Local Arrays ..
      DOUBLEPRECISION Cspl((NrDAT+1)*(NtDAT+1)),Fspl((NrDAT+1)*(NtDAT+
     *	1)),FF(NzDAT+1),LAMDA(NrDAT+5),MU(NtDAT+5),WRK(LWRK),Xrad(0:
     *    NzDAT),Ythet(0:NzDAT)
      INTEGER IWRK(LIWRK)
C-     .. External Subroutines ..
      EXTERNAL E01DAF,E02DEF
C-     .. Intrinsic Functions ..
      INTRINSIC MAX,MIN
      Integer mmax,Imax0,ii,im,ip,NtableIB,nDersMAX,nDersMX,nDers1,Narg,
     *	KMax,nnzero,LICN,LIRN,NOUT,L,irho0,IY,N_max,nStep,Mcheb,ITMIN
     *    ,ITMAX,IRES,krho,itable
      doubleprecision Xmin,Xmax,Zcyl,CMGM,expand,x_min,x_max,ploty_min,
     *	ploty_min2,ploty_max,ploty_max2,plotx,ploty,ploty2
      parameter(mmax=1000)
      dimension plotx(0:mmax),ploty(0:mmax),ploty2(0:mmax)
C- for sm
      character*80 filein,psfile,device,titlestr,xlabel,ylabel,string,
     *	charg
      Parameter(IMax0=1000)
      doubleprecision EOS,rcylm,rcylp,xm,xp,C,C0,c1,CM,V,dV,CJ,CI,T,W,
     *	Eth,Pip,VT,rhomax,H,P,psi,psiOut,omegaA
C-        double precision r, dr, dr0, Re, alpha, dalpha, theta, dtheta, phi, phi0, phiu, phiA,
C-          epsMain, rho, vphi, vphi2, Fdpsi_drcyl;
C-        Dimension r(0:IMax),alpha(0:IMax),theta(0:JMax);
C-      real etime,dtime,tarray(2) ! Convex
C-        Common /phi/phi(0:IMax,0:JMax);
C-        double precision rho0,rho01,rho02,delta,gamma,CN,CK;
C-        Common /delta/delta;
C-        Common /gamma/gamma,CN,CK;
C-        Common /pi/pi;
C-        Common /IB/IB/JB/JB/I0/I0;
      doubleprecision phi0
      Dimension phi0(0:IMax0)
      Parameter(NtableIB=10)
      Real*8 tableIB(NtableIB)
      data tableIB/1.,.917,.833,.750,.667,.625,.620,.619,.618,.617/
C- Hachisu
C-       Real*8 tableIB(NtableIB)/1.,.9,.8,.75,.7,.65,.64,.635,.63,.625/; -- AA290(1994)674
      Logical LogIter,LogInit,LogRapid,Lphi,Lrho,Lpolytrope
      Character*80 str,choblatness
      Character*6 rotlaw
C- rigid, jmscf, vconst, jconst or alpha
      REAL*4 secnds,t0,t1
C- timer
C-       LogRapid=.true.; -- default for CPsi in new_phi -- alpha-law rotation
      LogRapid=.false.
C- must be for other rotation laws test
C:  read rotlaw to be used  *
C------- '-->Entering Node %_arguments:'
      write(*,*)' Arguments needed: rotlaw -- either rigid,  
     *	jmscf, vconst, jconst, alpha or collap'
      write(*,*)' 2nd argument: oblatness: 0.< oblatness <= 1. '
      Narg=Iargc()
      GOTO(09999,09998,09997),Narg+1
      write(*,*)' Extra arguments.'
      Stop 16
      GOTO 09996
09999 CONTINUE ! Here the arguments are fixed
      rotlaw='rigid'
      oblatness=0.36d0  ! OBLATNESS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      GOTO 09996
09998 CONTINUE
      call GetArg(1,rotlaw)
      oblatness=0.8d0
      GOTO 09996
09997 CONTINUE
      call GetArg(1,rotlaw)
      oblatness=1.d0
      call GetArg(2,choblatness)
      read(choblatness,*)oblatness
09996 CONTINUE
      write(*,'(3a)')'  rotlaw=',rotlaw(1:len_trim(rotlaw)),' OK?'
      write(*,'(a,1p,1g11.3,a)')' oblatness=',oblatness,' OK?'
      call mypause
      if(rotlaw(1:len_trim(rotlaw)).EQ.'alpha')then
      LogRapid=.true.
      else
      LogRapid=.false.
C- LogRapid is always true when rotlaw is not alpha
      endif
C------- '<--Leaving  Node %_arguments:'
C- control vars, better read them from a dat file:
      JB=0
C- spheroidal configurations
      Lrho=.true.
C- the set of configurations along rho0 sequence
C-       Lrho=.false.; -- the set of configurations along IB  sequence
      Lphi=.true.
C- initiate phi
C-      Lphi=.false.; -- read phi from file phi
C: constants and (r,theta) grid *
C------- '-->Entering Node %_init:'
      NOUT=16
      open(1,file='rho0yudinm80.out',status='unknown')
      open(2,file='f.dat')
      open(3,file='inputSeb.dat')
      open(4,file='rho0MJ.cur',access='append')
      open(7,file='rCJMccyl.out',status='unknown')
C- this is for rho centered
      open(8,file='rCJMcylgrid.out',status='unknown')
C- this is for rho on cyl.grid
      open(NOUT,file='rho0Cheb.out',status='unknown')
      open(12,file='Prseb.out',status='unknown')
      read(3,'(f15.10,a80)')delta
C- grid parameter 1 or 2, true radius is r**delta,
C- Aksenov used delta=2
      read(3,'(f15.10,a80)')gamma
C-=4d0/3d0;
C-       gamma=5d0/3d0;
C-        gamma=2d0;
      CN=one/(gamma-one)
      CK=one/(one+CN)
C-       rho0=2d14; -- for neutron stars
      read(3,'(e20.10,a80)')rho0
C-=2.d09; -- for polytropic and eosseb
      read(3,'(L20,a80)')Lpolytrope
C- true for polytropic case
      close(3)
      LinitMcyl=.true.
C- mcyl must be initiated for jmscf rotating law
      write(*,'(a,1p,3g15.8)')'delta,gamma,rho0',delta,gamma,rho0
C-        pause;
C-      Call timer(ITime0) ! PC
      t0=secnds(0.0)
      write(1,'(18a)')'      rho0          ','r(IB)**delta','     Re
     *	         ','       CM_u/S_mass ','     CJ_u/1d50 ','    
     * T/abs(W) ','       gamma_  ',	'       rhomax  ','      phi0   '
     *,'        phiA   ','        phiB    ','        CM     ','        
     *V     ','        CJ     ','        CI    ','         W     ',
     *'         VT     ','        IB   '
C-                       '  trun';
      write(2,'(3a)')'         rho0,           r(IB)**delta,      Re,','
     *	             Mass,             V,           CJ,             
     *  CI,','       T/abs(W),       Pip/Eth+1,         VT'
      write(12,'(7a)')'#    r(i)   ','  rho(i,0)  ','     P      ',
     *	'     H      ','  phi(i,0)  ','    psi     ','  vphi'
      dr=one/IA
      do i=0,IA
      r(i)=dr*i
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
      enddo
      dalpha=one/r(IA-1)-one/r(IA)
      do l=0,IMax-IA+1
      alpha(l)=dalpha*l
      enddo
      do i=IA+1,IMax-1
      l=IMax-i
      r(i)=one/alpha(l)
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
      enddo
      dtheta=half*pi/dble(JMax)
      do j=0,JMax
      	theta(j)=dtheta*j
      enddo
      do ir=0,NrDAT
      Rcyld(ir)=powrsmart(r(ir),delta)
      enddo
C------- '<--Leaving  Node %_init:'
C:   initial phi - potential  *
C------- '-->Entering Node %_initphi:'
C-     === initial phi ==>
      if(Lphi)then
      do j=0,JMax
      do i=0,IMax
      if(i.le.IA)then
      phi(i,j)=float(i-IA-1)/float(IA)
      else
      l=IMax-i
      phi(i,j)=phi(IA,j)*powrsmart(alpha(l),delta)
      endif
      enddo
      enddo
      else
      open(16,file='phi',status='old')
      read(16,*)phi0
      close(16)
      do j=0,JMax
      do i=0,IMax
      if(i.le.IA)then
      dr0=one/IMax0
      im=powrsmart(r(i),delta)/dr0
      ip=im+1
      rcylm=im*dr0
      rcylp=ip*dr0
      xm=(rcylp-powrsmart(r(i),delta))/dr0
      xp=(powrsmart(r(i),delta)-rcylm)/dr0
      phi(i,j)=xm*phi0(im)+xp*phi0(ip)
      else
      l=IMax-i
      phi(i,j)=phi(IA,j)*powrsmart(alpha(l),delta)
      endif
      enddo
      enddo
      endif
C-     <== initial phi ===
C------- '<--Leaving  Node %_initphi:'
Cinitmcyl : initial mcyl()  *
      if(JB.eq.0)then
C- spheroidal configurations
      I0=0
      else
C- toroidal configurations
      I0=half*(IB+IA)
      endif
      if(Lrho)then
      rho01=1.D14   !   DEN1----------------------------------------------------------------------
      rho02=1.D13   !   DEN2----------------------------------------------------------------------
C- NS -- add to arguments
C-            rho01=5d4;
C           rho01=1d9;
C           rho02=5d9;  -- WD*
C-            irho0max=60; -- 20; -- not more than 999 for correct names of plot files
C-            write(4,'(a)')'  rho//J            M';
      do irho0=0,irho0max
C- Aksenov: ===== rho0 cycle ====>
C-                IB=IA; -- if IB==IA we have spherical configurations
C- IB=IA-80; -- if IB==IA we have spherical configurations
C- how to make non-zero rotation?
C-               IB=0.6*IA; -- try rotating sequence along rho0
      IB=oblatness*IA
C- assuming irho0<=9999
      IF(.NOT.(irho0.LT.10))GOTO 09995
      write(str,'(a,I1)')'000',irho0
      GOTO 09992
09995 CONTINUE
      IF(.NOT.(irho0.LT.100))GOTO 09994
      write(str,'(a,I2)')'00',irho0
      GOTO 09992
09994 CONTINUE
      IF(.NOT.(irho0.LT.1000))GOTO 09993
      write(str,'(a,I3)')'0',irho0
      GOTO 09992
09993 CONTINUE
      write(str,'(I4)')irho0
09992 CONTINUE
      open(26+irho0,file='phirho'//str(1:len_trim(str))//'.res')
      rho0=rho01+dble(irho0)/max(one,dble(irho0max))*(rho02-rho01)
      alphaPolytrope=sqrt((one+CN)*CK*rho0**(1.d0/CN-1.d0)/
     *	(4.d0*pi*G_N))
      write(*,'(a,1p,g12.3)')' rho0=',rho0
C-                call mypause;
C: Iterate Phi, C  *
C------- '-->Entering Node %_FindPhi:'
      write(*,*)' entering _FindPhi'
C-            call mypause;
      LogInit=.true.
      write(*,'(a,1p,g15.5,i7)')'rho0, IB',rho0,IB
      write(*,'(a,1p,3g15.5)')'phi(IA,JMax),phi(I0,JMax) =',phi(IA,JMax)
     *	,phi(I0,JMax)
C-           'phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)))=',phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)));
C- write of Fpsi is forbidden if there is write in it
      write(*,'(2a)')' rotlaw=',rotlaw(1:len_trim(rotlaw))
      C=(phi(IA,JMax)+Fpsi(one,rotlaw(1:len_trim(rotlaw)))-phi(I0,JMax)*
     *	EOS(zero,3)/EOS(rho0,3)
C- EOS(*,3) means enthalpy H
     *)/(one-EOS(zero,3)/EOS(rho0,3))
C- eq.17 in AA290(1994)674
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      write(*,'(a,1p,3g15.5)')'phiu,EOS(rho0,3),C=',phiu,EOS(rho0,3),C
C-            call mypause;
      n_max=300
      epsMain=1.d-12
C- to data file
      LogIter=.true.
      nStep=1
      do while(LogIter.and.(nStep.LE.n_max))
      C0=C
      Call new_phi(
C-r,alpha,theta,dtheta, -- now in module cylinderGrid
     *C,LogInit,LogRapid,rotlaw(1:len_trim(rotlaw)))
C- last argument is the rotation law
C              if(C>C0)then;
C                C=C0+min(abs(C-C0),0.03d0*abs(C0));
C              else;
C                C=C0-min(abs(C-C0),0.03d0*abs(C0));
C              endif;*
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      if(abs((C0-C)/C).LT.epsMain)LogIter=.false.
      write(*,*)' nStep,C,phi(I0,JMax):',nStep,C,phi(I0,JMax)
      nStep=nStep+1
      enddo
C------- '<--Leaving  Node %_FindPhi:'
C: to channel 26+irho0 after convergence, find integrals,
C                                remap onto cylindrical coordinates,
C                                prepare plots etc. *
C------- '-->Entering Node %_OutputResults:'
      Call Integrals(CM,V,CJ,CI,T,W,Eth,Pip,C,omegaA,
     *	rotlaw(1:len_trim(rotlaw)))
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' VT=',VT
      rhomax=zero
      do i=0,IA
C- loop on r
      do j=0,JMax
C- loop on theta
      if(j.EQ.JMax)then
      phiEqu(i)=phi(i,j)
C                   write(*,'(a,2i5,1p,2g15.6)')' i, j, theta(j),phiEqu(i)=',
C                            i,j,theta(j),phiEqu(i);
C                   if(mod(i,20)==0)call mypause;*
      endif
C- write(*,*)i,j,delta,C;
      rho=EOS(phiu*(C-phi(i,j)-Fpsi(powrsmart(r(i),delta)*sin(theta(j))
     *	,rotlaw(1:len_trim(rotlaw)))),1)
C- in phys. units from EOS
      rhoOn_r_thetaGrid(i,j)=rho
C- save for future
      rhomax=max(rhomax,rho/rho0)
      if(i.GT.0.and.j.EQ.Jmax.and.irho0.EQ.irho0max)then
      H=EOS(rho,3)
      P=EOS(H,4)
      psi=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))
      vphi2=-Fdpsi_drcyl(powrsmart(r(i),delta)*sin(theta(j)),
     *rotlaw(1:len_trim(rotlaw)))*powrsmart(r(i),delta)*sin(theta(JMax))
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,vphi2))
      write(12,'(1p,7g12.4)')r(i),rho,P,H,phi(i,j),psi,vphi
      endif
C: write to channel 26+irho0 for plots and debug
C                              in spherical coordinates  *
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
      psiOut=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))
C-                  Phirthet(i,j) = C -phi(i,j)-psiOut;
      Phirthet(i,j)=phi(i,j)
      write(26+irho0,'(1p,3e20.12)')powrsmart(r(i),delta),theta(j),
     *	C-phi(i,j)-psiOut
C
C  if( (i==IA .or. i==IA/2) .and. (j==1 .or. j==JMax) ) then; -- for debug
C        write(*,'(a,3i5,1p4g12.3)')' irho0 i j C phi Fpsi H:',irho0,
C                            i,j,C, phi(i,j), psiOut,C -phi(i,j)-psiOut;
C  write(*,'(a,1p,4e20.12)')' on write delta, r^delta, theta :',
C                    delta, powrsmart(r(i), delta), theta(j),
C                     C -phi(i,j)-psiOut;
C       call mypause;
C  endif;
C*
Coutcylinder : write to channel 26+IB for plots and debug
C                              in cylinder coordinates  *
      enddo
      enddo
C-            Phirthet(0,0) = C -phi(0,0);
      write(*,*)'Phirthet 00 11:',Phirthet(0,0),Phirthet(1,1)
C-            call mypause;
C-          open(1,file='output',access='append')
C-          write(1,*) 'time',etime(tarray),dtime(tarray)
C- === PC ==>
C-          Call timer(ITime)
C-           IDTime=ITime-ITime0;
      t1=secnds(0.0)-t0
C- <== PC ===
C-           write(1,*) '*** IB, time', IB, IDTime;
C-            write(1,'(a,i6,1p,e11.3)') '#*** IB, time:', IB, t1;
      Re=sqrt(phiu/(G_N*rho0))
      if(Lrho)then
      rhokeep(irho0)=rho0
      Jkeep(irho0)=CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50
      Mkeep(irho0)=CM*Re**3*rho0/S_mass
      endif
      write(*,*)'Re=',Re
C-            call mypause;
      write(1,22)rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,CJ*sqrt(G_N)
     *	*Re**5*rho0**1.5d0/1d50,T/abs(W),Pip/Eth+one,rhomax,phi(I0,JM
     *ax),phi(IA,JMax),phi(IB,JB),CM,V,CJ,CI,W,VT,IB
C-,t1;
C-          close(1)
      write(*,*)' output integrals'
      write(2,2)'intg',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     * ,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
C-            close(2);
      write(*,2)'intg',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
      write(*,*)'CJ, VT',CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,VT
C-            pause;
C: prepare spline for changing to cylindrical coords *
C------- '-->Entering Node %_OutputResults_Prepspline:'
C-        NrDAT=NrDAT;
C-        NzDAT=NzDAT;
      write(*,*)' N=',NrDAT,'   NzDAT=',NzDAT,'   NtDAT=',NtDAT
      write(*,*)' NrDAT=',NrDAT,'   NzDAT=',NzDAT
C-        call mypause;
      do ir=0,NrDAT
      rspl(ir)=powrsmart(r(ir),delta)
      enddo
      do jt=0,NtDAT
      thet(jt)=min(theta(jt),pi/2.d0)
      enddo
C-     .. Executable Statements ..
      WRITE(*,*)'E01DAF Program Results'
C-     the number of X points, MX, and the values of the
C-     X co-ordinates.
C-        MX=NrDAT;
      MX=NrDAT+1
C-     the number of Y points, MY, and the values of the
C-     Y co-ordinates.
      MY=NtDAT+1
C-     Read the function values at the grid points.
C       do jt=0,MY;
C          do ir=0,MX;
C            Fspl((MY+1)*ir+jt+1)=Phirthet(ir,jt);
C          enddo;
C       enddo;*
      do jt=1,MY
      do ir=1,MX
      Fspl(MY*(ir-1)+jt)=Phirthet(ir-1,jt-1)
      enddo
      enddo
      IFAIL=0
C-
C-  *     Generate the (X,Y,F) interpolating bicubic B-spline.
C-       CALL E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
C-     Generate the (rspl,thet,Fspl) interpolating bicubic B-spline.
      CALL E01DAF(MX,MY,rspl,thet,Fspl,PX,PY,LAMDA,MU,Cspl,WRK,IFAIL)
C- i.e. here the function Fspl(rspl,thet) is given on a grid NrDAT+1 \times NtDAT+1
C- as an array, the knots are in LAMDA and MU 1D arrays:
C
C      SUBROUTINE E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
CC     MARK 14 RELEASE. NAG COPYRIGHT 1989.
CC     Derived from DASL routine B2IRE.
CC     .. Parameters ..
C      CHARACTER*6       SRNAME
C      PARAMETER         (SRNAME='E01DAF')
C      DOUBLE PRECISION  ONE
C      PARAMETER         (ONE=1.0D+0)
CC     .. Scalar Arguments ..
C      INTEGER        IFAIL, MX, MY,
C         PX, PY -- output
CF(MX*MY) – real array  -- Input
COn entry: F(my * (q - 1)+r) must contain f_{q,r}, for q = 1, 2, ... ; mx , r = 1, 2, ... , my.
C
CC(MX*MY) – real array -- Output
COn exit: the coefﬁcients of the spline interpolant. C(my * (i - 1) + j)
Ccontains the coefﬁcient c_{ij} described in Section 3.
C
C
CC     .. Array Arguments ..
C      DOUBLE PRECISION  C(MX*MY), F(MX*MY), LAMDA(MX+4), MU(MY+4),
C     *                  WRK((MX+6)*(MY+6)), X(MX), Y(MY)
C_________________________________________________________________________
C*     E01DAF Example Program Text
C*     Mark 14 Release.  NAG Copyright 1989.
C*     .. Parameters ..
C      INTEGER          NIN, NOUT
C      PARAMETER        (NIN=5,NOUT=6)
C      INTEGER          MXMAX, MYMAX
C      PARAMETER        (MXMAX=20,MYMAX=MXMAX)
C      INTEGER          LIWRK, LWRK
C      PARAMETER        (LIWRK=MXMAX+2*(MXMAX-3)*(MYMAX-3),LWRK=(MXMAX+6)
C     +                 *(MYMAX+6))
C*     .. Local Scalars ..
C      DOUBLE PRECISION STEP, XHI, XLO, YHI, YLO
C      INTEGER          I, IFAIL, J, MX, MY, NX, NY, PX, PY
C*     .. Local Arrays ..
C      DOUBLE PRECISION C(MXMAX*MYMAX), F(MXMAX*MYMAX), FG(MXMAX*MYMAX),
C     +                 LAMDA(MXMAX+4), MU(MYMAX+4), TX(MXMAX),
C     +                 TY(MYMAX), WRK(LWRK), X(MXMAX), Y(MYMAX)
C      INTEGER          IWRK(LIWRK)
C      CHARACTER*10     CLABS(MYMAX), RLABS(MXMAX)
C  *
C-
C --     Print the knot sets, LAMDA and MU.
C       WRITE (*,*)'               I   Knot LAMDA(I)      J     Knot MU(J)';
C       do  J = 4, MAX(PX,PY) - 3;
C          IF (J.LE.PX-3 .AND. J.LE.PY-3) THEN;
C             WRITE (*, 99997) J, LAMDA(J), J, MU(J);
C          ELSE IF (J.LE.PX-3) THEN;
C             WRITE (*, 99997) J, LAMDA(J);
C          ELSE IF (J.LE.PY-3) THEN;
C             WRITE (*, 99996) J, MU(J);
C          END IF;
C       enddo;*
C-        call mypause;
C-     Print the spline coefficients.
C-        WRITE (*,*) 'The B-Spline coefficients computed in Cspl';
C-        WRITE (*,99999) (Cspl(I),I=1,MX*MY);
C------- '<--Leaving  Node %_OutputResults_Prepspline:'
C:  FROM Phirthet interpolate to phiRcylZ on a grid of Rcylindrical and Z *
C------- '-->Entering Node %_OutputResults_RcylZ:'
C- X(1)=0.; -- is Rcyl coord
C now in %_init:
C   do icyl=1,NrDAT;
C     Rcyld(icyl)=rspl(icyl);
C--      X(i)=Rcyld(i); -- X(i) not needed
C   enddo;
C   Rcyld(0)=0.d0;
C*
C- Y(1)=0.; -- is Z coord
      Zd(0)=0.d0
      do jcyl=1,NzDAT
      Zd(jcyl)=(Rcyld(NrDAT)/real(NzDAT))*real(jcyl)
      write(*,*)' jcyl,  Zd(jcyl)=',jcyl,Zd(jcyl)
C-       Y(jcyl)=Zd(jcyl); -- Y(i) not needed
      enddo
      do icyl=0,NrDAT
      do jcyl=0,NzDAT
C- height loop
      radius=sqrt(Rcyld(icyl)**2+Zd(jcyl)**2)
      if(jcyl.GT.0.d0)then
      thetspl=atan(Rcyld(icyl)/Zd(jcyl))
      else
      thetspl=pi/2.d0
      endif
C-       write(*,'(a,2i4,a,1p,4e11.3)')' icyl jcyl=',icyl,jcyl,' X Y radius thetspl =',X(icyl),
C-          Y(jcyl),radius,thetspl;
Clin : find Phi for point icyl,jcyl linearly interpolating from
C                 Phirthet at radius,thetspl *
      Xrad(jcyl)=max(rspl(0),min(radius,rspl(NrDAT)))
      Ythet(jcyl)=max(thet(0),min(thetspl,thet(NtDAT)))
C- call mypause;
      enddo
C- jcyl, i.e. height loop
C-      Xrad(0)=0.d0;
C-      Ythet(0)=0.d0;
C: find PhiRcylZ(icyl,jcyl) for NrDAT points icyl, and current jcyl by spline interpolation from
C                 Phirthet at radius,thetspl *
C-
C-        Evaluate the spline at M points
      M=NzDAT+1
      CALL E02DEF(M,PX,PY,Xrad(0),Ythet(0),LAMDA,MU,Cspl,FF,WRK,IWRK,
     *	IFAIL)
C
C      SUBROUTINE E02DEF(M,PX,PY,X,Y,LAMDA,MU,Cspl,FF,WRK,IWRK,IFAIL)
CC     MARK 14 RELEASE. NAG COPYRIGHT 1989.
CC
CC     Derived from DASL routine B2VRE.
CC
CC     E02DEF. An algorithm for evaluating a bicubic polynomial
CC     spline S(X,Y) from its B-spline representation at the M
CC     points (X(I),Y(I)), I = 1, 2, ..., M.
CC
CC     Input Parameters:
CC        M        The number of evaluation points.
CC        PX       NXKNTS + 8,  where  NXKNTS  is number
CC                    of interior  X-knots.
CC        PY       NYKNTS + 8,  where  NYKNTS  is number
CC                    of interior  Y-knots.
CC        X        X-values.
CC        Y        Y-values.
CC        LAMDA    The X-knots.
CC        MU       The Y-knots.
CC        Cspl        B-spline coefficients of  S.  First
CC                    subscript relates to  X.
CC
CC     Output parameter:
CC        FF       Values of spline.
CC                 On exit, FF(I) contains the value of
CC                 the spline evaluated at point (X(I),Y(I)),
CC                 for I = 1,..,M.
CC
CC     Workspace parameters:
CC        WRK      Real workspace of dimension at least (PY-4).
CC        IWRK     Integer workspace of dimension at least (PY-4).
CC
CC     Failure indicator parameter:
CC        IFAIL    Failure indicator:
CC                 1 -  PX     .LT. 8,    or
CC                      PY     .LT. 8,    or
CC                      M      .LT. 1.
CC                 2 -  E02DFW  failure for  X  or  Y.
CC                 3 -  At least one point (X(K),Y(K)) lies
CC                      outside the rectangle defined by
CC                      LAMDA(4), LAMDA(PX-3), MU(4) and MU(PY-3).
C  *
      If(Ifail.ne.0)then
      write(*,*)' Ifail after spline E02DEF=',Ifail
      call mypause
      endif
      do jcyl=1,NzDAT+1
      phiRcylZ(icyl,jcyl-1)=FF(jcyl)
C- this is dp either phi or C -phi(icyl,jcyl)-psiOut on Rcyl, Z grid
      if(icyl.GT.1.and.jcyl.GT.1)phiRcylZplot(icyl-1,jcyl-1)=FF(jcyl)
C- this is real*4 phi or C -phi(icyl,jcyl)-psiOut on Rcyl, Z grid
      enddo
      enddo
C- icyl, i.e. rcyl loop
C------- '<--Leaving  Node %_OutputResults_RcylZ:'
C: plots and/or control integrals in cylindrical coords *
C------- '-->Entering Node %_OutputResults_controls:'
C -- for supermongo
C        open(3,file='PhiRcylZ.res',form='unformatted'); -- this file is for sm image
C--  now write your data into arr
C        write(3) NrDAT,NzDAT;
C        write(3) phiRcylZ;
C*
C- -- for pgplot - pgxtal
C-      CALL PLOT(X,NrDAT,Y,NzDAT,PhiRcylZ,NrDAT,NzDAT,Wdummy,1.5,4,'Rcyl-axis','Z-axis',
C-                'potential phi(Rcyl,Z)');
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
      dZcyl=Zd(2)-Zd(1)
C- for uniform Zcyl grid
      write(*,*)' Zd(0)=',Zd(0)
      write(*,*)' dZcyl,Zd(1)-Zd(0)=',dZcyl,Zd(1)-Zd(0)
C-        call mypause;
C: for Zcyl=Zd(0) compare rho and phi in spherical
C                 and cylinder coordinates in the same points*
C------- '-->Entering Node %_OutputResults_controls_equator:'
      write(*,*)' NrDAT,IA=',NrDAT,IA
C: try different 1D interpolation procedures from NAG to produce phiEqu(Rcyl)
C      for z=0 and plot them online with supermongo  *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate:'
      write(*,*)' before OutputResults_controls_equator_interpolate
     *	_Chebyshev'
C-       call mypause;
C: use E01AEF to build Chebyshev interpolant and E02AKF to find interpolated values
C                  with  Chebyshev polynomials
C                  1) along z=0 for intermediate points of Rcyl grid;
C                  2) use phi on r,theta grid for constant r and interpolate phi(theta)  to theta=pi/2  *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate_Chebyshev:'
C-      E01AEF. A ROUTINE, WITH CHECKS, WHICH DETERMINES AND
C-      REFINES A POLYNOMIAL INTERPOLANT  Q(Xch)  TO DATA WHICH
C-      MAY CONTAIN DERIVATIVES.
C-
C-      INPUT PARAMETERS
C-         Mcheb        NUMBER OF DISTINCT Xch-VALUES
C-         XMIN,
C-         XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C-         Xch        INDEPENDENT VARIABLE VALUES (DISTINCT)
C-         Ych        VALUES AND DERIVATIVES OF
C-                     DEPENDENT VARIABLE.
C-         nDers       HIGHEST ORDER OF DERIVATIVE AT EACH Xch-VALUE.
C-         in the Rotat code we have no derivatives, so all nDers==0
C-         Ncheb        NUMBER OF INTERPOLATING CONDITIONS.
C-                     Ncheb == Mcheb + nDers(1) + nDers(2) + ... + nDers(Mcheb).
C-         in the Rotat code we have no derivatives, so Ncheb == Mcheb
C-         ITMIN,
C-         ITMAX    MINIMUM AND MAXIMUM NUMBER OF ITERATIONS TO BE
C-                     PERFORMED.
C-
C-      OUTPUT PARAMETERS
C-         chebA        CHEBYSHEV COEFFICIENTS OF  Q(Xch)
C-
C-      WORKSPACE (AND ASSOCIATED DIMENSION) PARAMETERS
C-         ChWRK      REAL WORKSPACE ARRAY.  THE FIRST IMAX ELEMENTS
C-                     CONTAIN, ON EXIT, PERFORMANCE INDICES FOR
C-                     THE INTERPOLATING POLYNOMIAL, AND THE NEXT
C-                     Ncheb  ELEMENTS THE COMPUTED RESIDUALS
C-         ChLWRK     DIMENSION OF ChWRK. ChLWRK MUST BE AT LEAST
C-                     7*Ncheb + 5*IMAX + Mcheb + 2, WHERE
C-                     IMAX IS ONE MORE THAN THE LARGEST ELEMENT
C-                     OF THE ARRAY nDers.
C-         IWRK     INTEGER WORKSPACE ARRAY.  ON EXIT,  IWRK(1)
C-                     CONTAINS THE NUMBER OF ITERATIONS TAKEN
C-         ChLIWRK    DIMENSION OF IWRK.  AT LEAST 2*Mcheb + 2.
C-
C-      FAILURE INDICATOR PARAMETER
C-         IFAIL    FAILURE INDICATOR.
C-                     0 - SUCCESSFUL TERMINATION.
C-                     1 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C-                         HAS BEEN VIOLATED -
C-                            Mcheb AT LEAST 1,
C-                            Ncheb = Mcheb + nDers(1) + nDers(2) + ... + nDers(Mcheb),
C-                            ChLWRK AT LEAST 7*Ncheb + 5*IMAX + Mcheb + 2,
C-                            ChLIWRK AT LEAST 2*Mcheb + 2.
C-                     2 - FOR SOME I, nDers(I) IS LESS THAN 0.
C-                     3 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C-                         HAS BEEN VIOLATED -
C-                            XMIN STRICTLY LESS THAN XMAX,
C-                            FOR EACH I, Xch(I) MUST LIE IN THE
C-                               INTERVAL XMIN TO XMAX,
C-                            THE Xch-VALUES MUST ALL BE DISTINCT
C-                     4 - NOT ALL PERFORMANCE INDICES LESS THAN
C-                         ONE, BUT ITMAX ITERATIONS PERFORMED,
C-                     5 - COMPUTATION TERMINATED BECAUSE
C-                         ITERATIONS DIVERGING.
C-
C-
C-      CHECK AND SET ITERATION LIMITS
C-
C-      .. Parameters ..
C-      CHARACTER*6       SRNAME
C-      PARAMETER         (SRNAME='E01AEF')
C-      .. Scalar Arguments ..
C-      DOUBLE PRECISION  XMAX, XMIN
C-      INTEGER           IFAIL, ITMAX, ITMIN, ChLIWRK, ChLWRK, Mcheb, Ncheb
C-      .. Array Arguments ..
C-      DOUBLE PRECISION  chebA(Ncheb), WRK(ChLWRK), Xch(Mcheb), Ych(Ncheb)
C-      INTEGER           nDers(Mcheb), IWRK(ChLIWRK)
      write(*,*)' entering OutputResults_controls_equator_
     *	interpolate_Chebyshev'
C: define for each I from 1 to Mcheb the values of nDers(I), Xch(I),
C                and (Ych(J),J=N+1,N+nDers(I)+1).
C            nDersMAX = MAX(nDersMAX,nDers(I))
C            N = N + nDers(I) + 1
C   *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate_Chebyshev_parChebyshev:'
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
      phia=phiRcylZ(icyl,0)
      enddo
      Mcheb=NrDAT+1
C- too many points
      Mcheb=8
C- try to do in pieces (may be overlapping), a loop over  pieces will be needed
      skipCheb=20
      Mcheb=(NrDAT+1)/skipCheb
      write(*,*)'  NrDAT+1,skipCheb,Mcheb=',NrDAT+1,skipCheb,Mcheb
C-          pause;
      Ncheb=0
      do I=1,Mcheb
      nDers(I)=0
C- non-zero if derivatives are given
      Xch(I)=Rcyld(skipCheb*I-1)
      do J=Ncheb+1,Ncheb+nDers(I)+1
      Ych(J)=phiRcylZ(skipCheb*J-1,0)
      enddo
      Ncheb=Ncheb+nDers(I)+1
      enddo
      XMIN=Xch(1)
      XMAX=Xch(Mcheb)
      IF(Ncheb.LE.ChNMAX.AND.nDersMAX.LE.nDerMX)THEN
      IFAIL=1
      nDersMAX=MAXVAL(nDers)
C- here must be zero
      write(*,*)' nDersMAX=',nDersMAX
C-             call mypause;
      ENDIF
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate_Chebyshev_parChebyshev:'
      write(*,*)' before   IF (Ncheb.LE.ChNMAX  ...)'
      write(*,*)' Ncheb.LE.ChNMAX .AND. nDersMAX.LE.nDersMX',Ncheb,
     *	ChNMAX,nDersMAX,nDersMX
C-          pause;
      IF(Ncheb.LE.ChNMAX.AND.nDersMAX.LE.nDersMX)THEN
      IFAIL=1
      ITMIN=0
      ITMAX=0
C- The next call to E01AEF can be used for obtaining  interpolant q(x) as a series of Chebyshev
C- polynomials with coefficients chebA(i).
C- The polynomial interpolant can subsequently be evaluated for any value of x in the
C- given range by using
C- E02AKF. Chebyshev-series representations of the derivative(s) and integral(s) of q(x) may be
C- obtained by (repeated) use of E02AHF and E02AJF.
      CALL E01AEF(Mcheb,XMIN,XMAX,Xch,Ych,nDers,Ncheb,ITMIN,ITMAX,chebA
     *	,ChWRK,ChLWRK,IWRK,ChLIWRK,IFAIL)
      write(*,*)' after call E01AEF(Mcheb, ...)'
C-             pause;
      WRITE(NOUT,*)
      IF(IFAIL.EQ.0.OR.IFAIL.GE.4)THEN
      WRITE(NOUT,89999)'Total number of interpolating conditions ='
     *	,Ncheb
      WRITE(NOUT,*)
      WRITE(NOUT,*)'Interpolating polynomial'
      WRITE(NOUT,*)
      WRITE(NOUT,*)'   I    Chebyshev Coefficient chebA(I+1)'
      DO I=1,Ncheb
      WRITE(NOUT,89998)I-1,chebA(I)
      ENDDO
      WRITE(NOUT,*)
      WRITE(NOUT,*)'  Xch    R   Rth derivative    Residual'
      IY=0
      IRES=nDersMAX+1
      DO I=1,Mcheb
      nDers1=nDers(I)+1
      DO J=1,nDers1
      IY=IY+1
      IRES=IRES+1
      IF(J-1.NE.0)THEN
      WRITE(NOUT,89997)J-1,Ych(IY),ChWRK(IRES)
      ELSE
      WRITE(NOUT,89996)Xch(I),'   0',Ych(IY),ChWRK(IRES)
      ENDIF
      ENDDO
      ENDDO
      ELSE
      WRITE(NOUT,89995)'E01AEF exits with IFAIL =',IFAIL
      write(*,*)' Ifail after Chebyshev E01AEF=',Ifail
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,Ncheb
     *	,XMIN,XMAX
      write(*,'(a,1p,(5g12.3))')' Xch(I)=',(Xch(I),I=1,MCheb)
C-                call mypause;
C-                pause; -- better than mypause, because exits smoothly on non-enter
      ENDIF
      WRITE(NOUT,89995)'Chebyshev E01AEF exits with IFAIL =',IFAIL
      write(*,*)' Ifail after Chebyshev E01AEF=',Ifail
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
      write(*,'(a,1p,(5g12.3))')' Xch(I)=',(Xch(I),I=1,MCheb)
      do i=1,MCheb
      IFAIL=0
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
C- this computes a function froom Chebyshev expansion:
      CALL E02AKF(NCheb+1,XMIN,XMAX,ChebA,1,NCheb+1,Xch(I),
     *	 resCheb(I),IFAIL)
      WRITE(*,'(a,i5,1p,2g12.3)')' ifail Xch(i) testCheb =',ifail,
     *	Xch(i),resCheb(i)
C- corresponds to ptype 6 3 in interactive sm:
      pt(i)=63.d0
      enddo
C-                call mypause;
C-                pause; -- better than mypause, because exits smoothly on non-enter
      ENDIF
C-       CALL E01BAF(N,Xch,Ych,LAMDA,C,LCK,WRK,LWRK,IFAIL);
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate_Chebyshev:'
C: for similar tasks E01BAF, or try the ready 2D spline *
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate:'
      Zcyl=Zd(0)
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
      phia=phiRcylZ(icyl,0)
C-                rho=EOS(phiu*(phiRcylZ(icyl,jcyl)),1)/rho0;
      psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)))
      rho=EOS(phiu*(C-phia-psi),1)
C- in phys. units from EOS
C-            write(*,'(a,1p,4g15.5)')' Rcyl, r, phia, psi=',Rcyl, r(icyl), phia, psi;
C-            write(*,'(a,1p,g15.5)')' rho on equ=',rho;
C-            write(*,'(a,1p,g15.5)')' rho saved =',rhoOn_r_thetaGrid(icyl,JMax);
C-       call mypause;
      enddo
C------- '<--Leaving  Node %_OutputResults_controls_equator:'
C: output CMcyl, CJcyl etc for rho centered *
C------- '-->Entering Node %_OutputResults_controls_rhocentered:'
      CMcyl(0)=0.d0
C- mass within Rcyl(icyl)
      CJcyl(0)=0.d0
C- angular momentum within Rcyl(icyl)
      write(7,'(3a)')'#     Rcyla,              CMcyl(icyl),           
     *	sigma','             sigmaGM,              ERROR','       IFA
     *IL,    CJcyl(icyl)'
      do icyl=1,NrDAT
      Rcyl=Rcyld(icyl)
C- here dp
      if(icyl.GT.1)then
C- now this "if  " is not needed
      Rcylm=Rcyld(icyl-1)
      else
      Rcylm=0.d0
      endif
      Rcyla=half*(Rcyl+Rcylm)
C-            psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)));
      psi=Fpsi(Rcyla,rotlaw(1:len_trim(rotlaw)))
      c1=-Fdpsi_drcyl(Rcyla,rotlaw(1:len_trim(rotlaw)))*Rcyla
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,c1))
      CMcyl(icyl)=CMcyl(icyl-1)
      CJcyl(icyl)=CJcyl(icyl-1)
      sigma=0.d0
      do jcyl=1,NzDAT
      Zcyl=Zd(jcyl)
C- here dp
C
C               if (jcyl==1) then;
C                  write(*,*)' Zcyl=',Zcyl;
C                  call mypause;
C               endif;
C*
      phia=(phiRcylZ(icyl-1,jcyl-1)+phiRcylZ(icyl-1,jcyl)+
     *	phiRcylZ(icyl,jcyl-1)+phiRcylZ(icyl,jcyl))/four
C-                rho=EOS(phiu*(phiRcylZ(icyl,jcyl)),1)/rho0;
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      rhoRcylZ(jcyl,icyl)=rho
C- 1st index for z
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-             call mypause;
      dV=2.d0*(Rcyl**2-Rcylm**2)*dZcyl*pi
C- multiplied by 2 to account for 'southern' hemisphere
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
C-  	       rhodV(jcyl)=2.d0*pi*rho*(Rcyl**2-Rcylm**2);
      	CM=CM+rho*dV
      	sigma=sigma+2.d0*rho*dZcyl
C- multiplied by 2 to account for 'southern' hemisphere
      CMcyl(icyl)=CMcyl(icyl)+rho*dV
      CJcyl(icyl)=CJcyl(icyl)+rho*vphi*Rcyla*dV
      	CJ=CJ+rho*vphi*Rcyla*dV
      	CI=CI+rho*Rcyla**2*dV
C-                P=EOS(phiu*(phiRcylZ(icyl,jcyl)),4)/(rho0*phiu);
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
C-  	       W=W+rho*phiRcylZ(icyl,jcyl)/2d0*dV;
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      enddo
C- loop on z
      Call GillMiller(Zd(1),rhoRcylZ(1,icyl),NzDAT,ANS,ERROR,IFAIL)
      sigmaGM=2.d0*ANS
C- multiplied by 2 to account for 'southern' hemisphere
      sigma2pir(icyl)=2.d0*pi*sigmaGM*Rcyla
      write(7,'(1p,5g20.8,i5,g20.8)')Rcyla,CMcyl(icyl),sigma,sigmaGM,
     *	ERROR,IFAIL,CJcyl(icyl)
      enddo
C- loop on Rcyl
      Call GillMiller(Rcyld(1),sigma2pir(1),NrDAT,CMGM,ERROR,IFAIL)
      write(7,'(a,1p,3g20.8,i5)')'# CM CMGM ERROR IFAIL=',CM,CMGM,
     *	ERROR,IFAIL
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' output controls'
      write(2,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     *	V*Re**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),
     *    Pip/Eth+one,VT
      write(*,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     *	V*Re**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),
     *    Pip/Eth+one,VT
      if(Lpolytrope)then
      callpolytrope(CN)
      write(2,2)'poly',rho0,r(IB)**delta/alphaPolytrope,Re/
     *	alphaPolytrope,CM*(Re/alphaPolytrope)**3/(4.d0*pi),V*(Re/
     *    alphaPolytrope)*
     **3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI*Re**5/(alphaPolytrope**5
     *),T/abs(W),Pip/Eth+one,VT
      endif
C-        pause;
C------- '<--Leaving  Node %_OutputResults_controls_rhocentered:'
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
C: output CMcylg, CJcylg etc for rho on cylindrical grid  *
C------- '-->Entering Node %_OutputResults_controls_rhoOnCylGrid:'
      CMcylg(0)=0.d0
C- mass within Rcyl(icyl)
      CJcylg(0)=0.d0
C- angular momentum within Rcyl(icyl)
      write(8,'(3a)')'#     Rcyla,              CMcylg(icyl),        
     *   sigma','             sigmaGM,              ERROR','       IF
     *AIL,    CJcylg(icyl)'
      Rcylm=0.d0
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
C- here dp
      psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)))
      c1=-Fdpsi_drcyl(Rcyl,rotlaw(1:len_trim(rotlaw)))*Rcyla
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,c1))
      if(icyl.GT.0)then
      CMcylg(icyl)=CMcylg(icyl-1)
      CJcylg(icyl)=CJcylg(icyl-1)
      Rcylm=Rcyld(icyl-1)
      endif
      sigma=0.d0
      dojcyl=0,NzDAT
      Zcyl=Zd(jcyl)
C- here dp
C
C               if (jcyl==1) then;
C                  write(*,*)' Zcyl=',Zcyl;
C                  call mypause;
C               endif;
C*
      phia=phiRcylZ(icyl,jcyl)
C- here phia on cylindrical grid
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      rhoRcylZ(jcyl,icyl)=rho
C- 1st index for z
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-             call mypause;
      dRcyl=Rcyl-Rcylm
      if(jcyl.LT.NzDAT)then
      dV=2.d0*((Rcyl**2-Rcylm**2)+dRcyl**2)*dZcyl*pi
C- multiplied by 2 to account for 'southern' hemisphere
      sigma=sigma+2.d0*rho*dZcyl
C- multiplied by 2 to account for 'southern' hemisphere
      else
      dV=((Rcyl**2-Rcylm**2)+dRcyl**2)*dZcyl*pi
C- not multiplied by 2 since dZcyl is on border
      sigma=sigma+rho*dZcyl
C- not multiplied by 2
      endif
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
C-  	       rhodV(jcyl)=2.d0*pi*rho*(Rcyl**2-Rcylm**2);
      	CM=CM+rho*dV
      CMcylg(icyl)=CMcylg(icyl)+rho*dV
      CJcylg(icyl)=CJcylg(icyl)+rho*vphi*Rcyla*dV
      	CJ=CJ+rho*vphi*Rcyla*dV
      	CI=CI+rho*Rcyla**2*dV
C-                P=EOS(phiu*(phiRcylZ(icyl,jcyl)),4)/(rho0*phiu);
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
C-  	       W=W+rho*phiRcylZ(icyl,jcyl)/2d0*dV;
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      enddo
C- loop on z
      Call GillMiller(Zd,rhoRcylZ(0,icyl),NzDAT+1,ANS,ERROR,IFAIL)
      sigmaGM=2.d0*ANS
C- multiplied by 2 to account for 'southern' hemisphere
      sigma2pir(icyl)=2.d0*pi*sigmaGM*Rcyl
      write(8,'(1p,5g20.8,i5,g20.8)')Rcyl,CMcylg(icyl),sigma,sigmaGM,
     *	ERROR,IFAIL,CJcylg(icyl)
      enddo
C- loop on Rcyl
      Call GillMiller(Rcyld,sigma2pir,NrDAT+1,CMGM,ERROR,IFAIL)
      write(8,'(a,1p,3g20.8,i5)')'# CM CMGM ERROR IFAIL=',CM,CMGM,
     *	ERROR,IFAIL
      do icyl=0,NrDAT
      mcyl(icyl)=CMcylg(icyl)/CMcylg(NrDAT)
      enddo
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' output controls'
      write(8,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     * 	V*Re**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),
     *    Pip/Eth+one,VT
      write(*,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     *	V*Re**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),
     *    Pip/Eth+one,VT
      if(Lpolytrope)write(8,2)'poly',rho0,r(IB)**delta/alphaPolytrope,
     *	Re/alphaPolytrope,CM*(Re/alphaPolytrope)**3/(4.d0*pi),V*(Re/al
     *phaPolytrope)**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI*Re**5
     */(alphaPolytrope**5),
C- *rho0),
     *T/abs(W),Pip/Eth+one,VT
C-        pause;
      LinitMcyl=.false.
C- needed for jmscf -- 1st model is rigid, 2nd jmscf
C------- '<--Leaving  Node %_OutputResults_controls_rhoOnCylGrid:'
C------- '<--Leaving  Node %_OutputResults_controls:'
C------- '<--Leaving  Node %_OutputResults:'
      close(26+irho0)
      enddo
C-            write(4,@ap(1p,@irho2max g15.5)@ap)((rhokeep(krho),rhokeep(krho)),krho=0,irho0max);
C-            write(4,@ap(1p,@irho2max g15.5)@ap)((Jkeep(krho),Mkeep(krho)),krho=0,irho0max);
      write(4,' (1p,102g15.5)' )(Jkeep(krho),Mkeep(krho),krho=0,
     *	irho0max)
      elseif(rotlaw.EQ.'rigid')then
C: read IB from tableIB as in Hachisu or in our paper
C                       AA290(1994)674 *
C------- '-->Entering Node %_tableIB:'
      do itable=1,NtableIB
      IB=nint(IA*tableIB(itable))
C- assuming IB<=9999
      IF(.NOT.(IB.LT.10))GOTO 09991
      write(str,'(a,I1)')'000',IB
      GOTO 09988
09991 CONTINUE
      IF(.NOT.(IB.LT.100))GOTO 09990
      write(str,'(a,I2)')'00',IB
      GOTO 09988
09990 CONTINUE
      IF(.NOT.(IB.LT.1000))GOTO 09989
      write(str,'(a,I3)')'0',IB
      GOTO 09988
09989 CONTINUE
      write(str,'(I4)')IB
09988 CONTINUE
      open(26+IB,file='phi'//str(1:len_trim(str))//'.res')
C: iterate phi etc. for rigid rotation *
C------- '-->Entering Node %_FindPhi:'
      write(*,*)' entering _FindPhi'
C-            call mypause;
      LogInit=.true.
      write(*,'(a,1p,g15.5,i7)')'rho0, IB',rho0,IB
      write(*,'(a,1p,3g15.5)')'phi(IA,JMax),phi(I0,JMax) =',
     *	phi(IA,JMax),phi(I0,JMax)
C-           'phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)))=',phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)));
C- write of Fpsi is forbidden if there is write in it
      write(*,'(2a)')' rotlaw=',rotlaw(1:len_trim(rotlaw))
      C=(phi(IA,JMax)+Fpsi(one,rotlaw(1:len_trim(rotlaw)))-phi(I0,JMax)
     *	*EOS(zero,3)/EOS(rho0,3)
C- EOS(*,3) means enthalpy H
     *)/(one-EOS(zero,3)/EOS(rho0,3))
C- eq.17 in AA290(1994)674
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      write(*,'(a,1p,3g15.5)')'phiu,EOS(rho0,3),C=',phiu,EOS(rho0,3),C
C-            call mypause;
      n_max=300
      epsMain=1.d-12
C- to data file
      LogIter=.true.
      nStep=1
      do while(LogIter.and.(nStep.LE.n_max))
      C0=C
      Call new_phi(
C-r,alpha,theta,dtheta, -- now in module cylinderGrid
     *C,LogInit,LogRapid,rotlaw(1:len_trim(rotlaw)))
C- last argument is the rotation law
C              if(C>C0)then;
C                C=C0+min(abs(C-C0),0.03d0*abs(C0));
C              else;
C                C=C0-min(abs(C-C0),0.03d0*abs(C0));
C              endif;*
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      if(abs((C0-C)/C).LT.epsMain)LogIter=.false.
      write(*,*)' nStep,C,phi(I0,JMax):',nStep,C,phi(I0,JMax)
      nStep=nStep+1
      enddo
C------- '<--Leaving  Node %_FindPhi:'
C: to channel 26+irho0 after convergence, find integrals,
C                                remap onto cylindrical coordinates,
C                                prepare plots etc. *
C------- '-->Entering Node %_OutputResults:'
      Call Integrals(CM,V,CJ,CI,T,W,Eth,Pip,C,omegaA,
     *	rotlaw(1:len_trim(rotlaw)))
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' VT=',VT
      rhomax=zero
      do i=0,IA
C- loop on r
      do j=0,JMax
C- loop on theta
      if(j.EQ.JMax)then
      phiEqu(i)=phi(i,j)
C                   write(*,'(a,2i5,1p,2g15.6)')' i, j, theta(j),phiEqu(i)=',
C                            i,j,theta(j),phiEqu(i);
C                   if(mod(i,20)==0)call mypause;*
      endif
C- write(*,*)i,j,delta,C;
      rho=EOS(phiu*(C-phi(i,j)-Fpsi(powrsmart(r(i),delta)*sin(theta(j))
     *	,rotlaw(1:len_trim(rotlaw)))),1)
C- in phys. units from EOS
      rhoOn_r_thetaGrid(i,j)=rho
C- save for future
      rhomax=max(rhomax,rho/rho0)
      if(i.GT.0.and.j.EQ.Jmax.and.irho0.EQ.irho0max)then
      H=EOS(rho,3)
      P=EOS(H,4)
      psi=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))
      vphi2=-Fdpsi_drcyl(powrsmart(r(i),delta)*sin(theta(j)),
     *rotlaw(1:len_trim(rotlaw)))*powrsmart(r(i),delta)*sin(theta(JMax))
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,vphi2))
      write(12,'(1p,7g12.4)')r(i),rho,P,H,phi(i,j),psi,vphi
      endif
C: write to channel 26+irho0 for plots and debug
C                              in spherical coordinates  *
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
      psiOut=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))
C-                  Phirthet(i,j) = C -phi(i,j)-psiOut;
      Phirthet(i,j)=phi(i,j)
      write(26+irho0,'(1p,3e20.12)')powrsmart(r(i),delta),theta(j),
     *	C-phi(i,j)-psiOut
C
C  if( (i==IA .or. i==IA/2) .and. (j==1 .or. j==JMax) ) then; -- for debug
C        write(*,'(a,3i5,1p4g12.3)')' irho0 i j C phi Fpsi H:',irho0,
C                            i,j,C, phi(i,j), psiOut,C -phi(i,j)-psiOut;
C  write(*,'(a,1p,4e20.12)')' on write delta, r^delta, theta :',
C                    delta, powrsmart(r(i), delta), theta(j),
C                     C -phi(i,j)-psiOut;
C       call mypause;
C  endif;
C*
Coutcylinder : write to channel 26+IB for plots and debug
C                              in cylinder coordinates  *
      enddo
      enddo
C-            Phirthet(0,0) = C -phi(0,0);
      write(*,*)'Phirthet 00 11:',Phirthet(0,0),Phirthet(1,1)
C-            call mypause;
C-          open(1,file='output',access='append')
C-          write(1,*) 'time',etime(tarray),dtime(tarray)
C- === PC ==>
C-          Call timer(ITime)
C-           IDTime=ITime-ITime0;
      t1=secnds(0.0)-t0
C- <== PC ===
C-           write(1,*) '*** IB, time', IB, IDTime;
C-            write(1,'(a,i6,1p,e11.3)') '#*** IB, time:', IB, t1;
      Re=sqrt(phiu/(G_N*rho0))
      if(Lrho)then
      rhokeep(irho0)=rho0
      Jkeep(irho0)=CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50
      Mkeep(irho0)=CM*Re**3*rho0/S_mass
      endif
      write(*,*)'Re=',Re
C-            call mypause;
      write(1,22)rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,CJ*sqrt(G_N)
     *	*Re**5*rho0**1.5d0/1d50,T/abs(W),Pip/Eth+one,rhomax,phi(I0,JM
     *ax),phi(IA,JMax),phi(IB,JB),CM,V,CJ,CI,W,VT,IB
C-,t1;
C-          close(1)
      write(*,*)' output integrals'
      write(2,2)'intg',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
C-            close(2);
      write(*,2)'intg',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
      write(*,*)'CJ, VT',CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,VT
C-            pause;
C: prepare spline for changing to cylindrical coords *
C------- '-->Entering Node %_OutputResults_Prepspline:'
C-        NrDAT=NrDAT;
C-        NzDAT=NzDAT;
      write(*,*)' N=',NrDAT,'   NzDAT=',NzDAT,'   NtDAT=',NtDAT
      write(*,*)' NrDAT=',NrDAT,'   NzDAT=',NzDAT
C-        call mypause;
      do ir=0,NrDAT
      rspl(ir)=powrsmart(r(ir),delta)
      enddo
      do jt=0,NtDAT
      thet(jt)=min(theta(jt),pi/2.d0)
      enddo
C-     .. Executable Statements ..
      WRITE(*,*)'E01DAF Program Results'
C-     the number of X points, MX, and the values of the
C-     X co-ordinates.
C-        MX=NrDAT;
      MX=NrDAT+1
C-     the number of Y points, MY, and the values of the
C-     Y co-ordinates.
      MY=NtDAT+1
C-     Read the function values at the grid points.
C       do jt=0,MY;
C          do ir=0,MX;
C            Fspl((MY+1)*ir+jt+1)=Phirthet(ir,jt);
C          enddo;
C       enddo;*
      do jt=1,MY
      do ir=1,MX
      Fspl(MY*(ir-1)+jt)=Phirthet(ir-1,jt-1)
      enddo
      enddo
      IFAIL=0
C-
C-  *     Generate the (X,Y,F) interpolating bicubic B-spline.
C-       CALL E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
C-     Generate the (rspl,thet,Fspl) interpolating bicubic B-spline.
      CALL E01DAF(MX,MY,rspl,thet,Fspl,PX,PY,LAMDA,MU,Cspl,WRK,IFAIL)
C- i.e. here the function Fspl(rspl,thet) is given on a grid NrDAT+1 \times NtDAT+1
C- as an array, the knots are in LAMDA and MU 1D arrays:
C
C      SUBROUTINE E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
CC     MARK 14 RELEASE. NAG COPYRIGHT 1989.
CC     Derived from DASL routine B2IRE.
CC     .. Parameters ..
C      CHARACTER*6       SRNAME
C      PARAMETER         (SRNAME='E01DAF')
C      DOUBLE PRECISION  ONE
C      PARAMETER         (ONE=1.0D+0)
CC     .. Scalar Arguments ..
C      INTEGER        IFAIL, MX, MY,
C         PX, PY -- output
CF(MX*MY) – real array  -- Input
COn entry: F(my * (q - 1)+r) must contain f_{q,r}, for q = 1, 2, ... ; mx , r = 1, 2, ... , my.
C
CC(MX*MY) – real array -- Output
COn exit: the coefﬁcients of the spline interpolant. C(my * (i - 1) + j)
Ccontains the coefﬁcient c_{ij} described in Section 3.
C
C
CC     .. Array Arguments ..
C      DOUBLE PRECISION  C(MX*MY), F(MX*MY), LAMDA(MX+4), MU(MY+4),
C     *                  WRK((MX+6)*(MY+6)), X(MX), Y(MY)
C_________________________________________________________________________
C*     E01DAF Example Program Text
C*     Mark 14 Release.  NAG Copyright 1989.
C*     .. Parameters ..
C      INTEGER          NIN, NOUT
C      PARAMETER        (NIN=5,NOUT=6)
C      INTEGER          MXMAX, MYMAX
C      PARAMETER        (MXMAX=20,MYMAX=MXMAX)
C      INTEGER          LIWRK, LWRK
C      PARAMETER        (LIWRK=MXMAX+2*(MXMAX-3)*(MYMAX-3),LWRK=(MXMAX+6)
C     +                 *(MYMAX+6))
C*     .. Local Scalars ..
C      DOUBLE PRECISION STEP, XHI, XLO, YHI, YLO
C      INTEGER          I, IFAIL, J, MX, MY, NX, NY, PX, PY
C*     .. Local Arrays ..
C      DOUBLE PRECISION C(MXMAX*MYMAX), F(MXMAX*MYMAX), FG(MXMAX*MYMAX),
C     +                 LAMDA(MXMAX+4), MU(MYMAX+4), TX(MXMAX),
C     +                 TY(MYMAX), WRK(LWRK), X(MXMAX), Y(MYMAX)
C      INTEGER          IWRK(LIWRK)
C      CHARACTER*10     CLABS(MYMAX), RLABS(MXMAX)
C  *
C-
C --     Print the knot sets, LAMDA and MU.
C       WRITE (*,*)'               I   Knot LAMDA(I)      J     Knot MU(J)';
C       do  J = 4, MAX(PX,PY) - 3;
C          IF (J.LE.PX-3 .AND. J.LE.PY-3) THEN;
C             WRITE (*, 99997) J, LAMDA(J), J, MU(J);
C          ELSE IF (J.LE.PX-3) THEN;
C             WRITE (*, 99997) J, LAMDA(J);
C          ELSE IF (J.LE.PY-3) THEN;
C             WRITE (*, 99996) J, MU(J);
C          END IF;
C       enddo;*
C-        call mypause;
C-     Print the spline coefficients.
C-        WRITE (*,*) 'The B-Spline coefficients computed in Cspl';
C-        WRITE (*,99999) (Cspl(I),I=1,MX*MY);
C------- '<--Leaving  Node %_OutputResults_Prepspline:'
C:  FROM Phirthet interpolate to phiRcylZ on a grid of Rcylindrical and Z *
C------- '-->Entering Node %_OutputResults_RcylZ:'
C- X(1)=0.; -- is Rcyl coord
C now in %_init:
C   do icyl=1,NrDAT;
C     Rcyld(icyl)=rspl(icyl);
C--      X(i)=Rcyld(i); -- X(i) not needed
C   enddo;
C   Rcyld(0)=0.d0;
C*
C- Y(1)=0.; -- is Z coord
      Zd(0)=0.d0
      do jcyl=1,NzDAT
      Zd(jcyl)=(Rcyld(NrDAT)/real(NzDAT))*real(jcyl)
      write(*,*)' jcyl,  Zd(jcyl)=',jcyl,Zd(jcyl)
C-       Y(jcyl)=Zd(jcyl); -- Y(i) not needed
      enddo
      do icyl=0,NrDAT
      do jcyl=0,NzDAT
C- height loop
      radius=sqrt(Rcyld(icyl)**2+Zd(jcyl)**2)
      if(jcyl.GT.0.d0)then
      thetspl=atan(Rcyld(icyl)/Zd(jcyl))
      else
      thetspl=pi/2.d0
      endif
C-       write(*,'(a,2i4,a,1p,4e11.3)')' icyl jcyl=',icyl,jcyl,' X Y radius thetspl =',X(icyl),
C-          Y(jcyl),radius,thetspl;
Clin : find Phi for point icyl,jcyl linearly interpolating from
C                 Phirthet at radius,thetspl *
      Xrad(jcyl)=max(rspl(0),min(radius,rspl(NrDAT)))
      Ythet(jcyl)=max(thet(0),min(thetspl,thet(NtDAT)))
C- call mypause;
      enddo
C- jcyl, i.e. height loop
C-      Xrad(0)=0.d0;
C-      Ythet(0)=0.d0;
C: find PhiRcylZ(icyl,jcyl) for NrDAT points icyl, and current jcyl by spline interpolation from
C                 Phirthet at radius,thetspl *
C-
C-        Evaluate the spline at M points
      M=NzDAT+1
      CALL E02DEF(M,PX,PY,Xrad(0),Ythet(0),LAMDA,MU,Cspl,FF,WRK,IWRK,
     *	IFAIL)
C
C      SUBROUTINE E02DEF(M,PX,PY,X,Y,LAMDA,MU,Cspl,FF,WRK,IWRK,IFAIL)
CC     MARK 14 RELEASE. NAG COPYRIGHT 1989.
CC
CC     Derived from DASL routine B2VRE.
CC
CC     E02DEF. An algorithm for evaluating a bicubic polynomial
CC     spline S(X,Y) from its B-spline representation at the M
CC     points (X(I),Y(I)), I = 1, 2, ..., M.
CC
CC     Input Parameters:
CC        M        The number of evaluation points.
CC        PX       NXKNTS + 8,  where  NXKNTS  is number
CC                    of interior  X-knots.
CC        PY       NYKNTS + 8,  where  NYKNTS  is number
CC                    of interior  Y-knots.
CC        X        X-values.
CC        Y        Y-values.
CC        LAMDA    The X-knots.
CC        MU       The Y-knots.
CC        Cspl        B-spline coefficients of  S.  First
CC                    subscript relates to  X.
CC
CC     Output parameter:
CC        FF       Values of spline.
CC                 On exit, FF(I) contains the value of
CC                 the spline evaluated at point (X(I),Y(I)),
CC                 for I = 1,..,M.
CC
CC     Workspace parameters:
CC        WRK      Real workspace of dimension at least (PY-4).
CC        IWRK     Integer workspace of dimension at least (PY-4).
CC
CC     Failure indicator parameter:
CC        IFAIL    Failure indicator:
CC                 1 -  PX     .LT. 8,    or
CC                      PY     .LT. 8,    or
CC                      M      .LT. 1.
CC                 2 -  E02DFW  failure for  X  or  Y.
CC                 3 -  At least one point (X(K),Y(K)) lies
CC                      outside the rectangle defined by
CC                      LAMDA(4), LAMDA(PX-3), MU(4) and MU(PY-3).
C  *
      If(Ifail.ne.0)then
      write(*,*)' Ifail after spline E02DEF=',Ifail
      call mypause
      endif
      do jcyl=1,NzDAT+1
      phiRcylZ(icyl,jcyl-1)=FF(jcyl)
C- this is dp either phi or C -phi(icyl,jcyl)-psiOut on Rcyl, Z grid
      if(icyl.GT.1.and.jcyl.GT.1)phiRcylZplot(icyl-1,jcyl-1)=FF(jcyl)
C- this is real*4 phi or C -phi(icyl,jcyl)-psiOut on Rcyl, Z grid
      enddo
      enddo
C- icyl, i.e. rcyl loop
C------- '<--Leaving  Node %_OutputResults_RcylZ:'
C: plots and/or control integrals in cylindrical coords *
C------- '-->Entering Node %_OutputResults_controls:'
C -- for supermongo
C        open(3,file='PhiRcylZ.res',form='unformatted'); -- this file is for sm image
C--  now write your data into arr
C        write(3) NrDAT,NzDAT;
C        write(3) phiRcylZ;
C*
C- -- for pgplot - pgxtal
C-      CALL PLOT(X,NrDAT,Y,NzDAT,PhiRcylZ,NrDAT,NzDAT,Wdummy,1.5,4,'Rcyl-axis','Z-axis',
C-                'potential phi(Rcyl,Z)');
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
      dZcyl=Zd(2)-Zd(1)
C- for uniform Zcyl grid
      write(*,*)' Zd(0)=',Zd(0)
      write(*,*)' dZcyl,Zd(1)-Zd(0)=',dZcyl,Zd(1)-Zd(0)
C-        call mypause;
C: for Zcyl=Zd(0) compare rho and phi in spherical
C                 and cylinder coordinates in the same points*
C------- '-->Entering Node %_OutputResults_controls_equator:'
      write(*,*)' NrDAT,IA=',NrDAT,IA
C: try different 1D interpolation procedures from NAG to produce phiEqu(Rcyl)
C      for z=0 and plot them online with supermongo  *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate:'
      write(*,*)' before OutputResults_controls_equator_interpolate
     *	_Chebyshev'
C-       call mypause;
C: use E01AEF to build Chebyshev interpolant and E02AKF to find interpolated values
C                  with  Chebyshev polynomials
C                  1) along z=0 for intermediate points of Rcyl grid;
C                  2) use phi on r,theta grid for constant r and interpolate phi(theta)  to theta=pi/2  *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate_Chebyshev:'
C-      E01AEF. A ROUTINE, WITH CHECKS, WHICH DETERMINES AND
C-      REFINES A POLYNOMIAL INTERPOLANT  Q(Xch)  TO DATA WHICH
C-      MAY CONTAIN DERIVATIVES.
C-
C-      INPUT PARAMETERS
C-         Mcheb        NUMBER OF DISTINCT Xch-VALUES
C-         XMIN,
C-         XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C-         Xch        INDEPENDENT VARIABLE VALUES (DISTINCT)
C-         Ych        VALUES AND DERIVATIVES OF
C-                     DEPENDENT VARIABLE.
C-         nDers       HIGHEST ORDER OF DERIVATIVE AT EACH Xch-VALUE.
C-         in the Rotat code we have no derivatives, so all nDers==0
C-         Ncheb        NUMBER OF INTERPOLATING CONDITIONS.
C-                     Ncheb == Mcheb + nDers(1) + nDers(2) + ... + nDers(Mcheb).
C-         in the Rotat code we have no derivatives, so Ncheb == Mcheb
C-         ITMIN,
C-         ITMAX    MINIMUM AND MAXIMUM NUMBER OF ITERATIONS TO BE
C-                     PERFORMED.
C-
C-      OUTPUT PARAMETERS
C-         chebA        CHEBYSHEV COEFFICIENTS OF  Q(Xch)
C-
C-      WORKSPACE (AND ASSOCIATED DIMENSION) PARAMETERS
C-         ChWRK      REAL WORKSPACE ARRAY.  THE FIRST IMAX ELEMENTS
C-                     CONTAIN, ON EXIT, PERFORMANCE INDICES FOR
C-                     THE INTERPOLATING POLYNOMIAL, AND THE NEXT
C-                     Ncheb  ELEMENTS THE COMPUTED RESIDUALS
C-         ChLWRK     DIMENSION OF ChWRK. ChLWRK MUST BE AT LEAST
C-                     7*Ncheb + 5*IMAX + Mcheb + 2, WHERE
C-                     IMAX IS ONE MORE THAN THE LARGEST ELEMENT
C-                     OF THE ARRAY nDers.
C-         IWRK     INTEGER WORKSPACE ARRAY.  ON EXIT,  IWRK(1)
C-                     CONTAINS THE NUMBER OF ITERATIONS TAKEN
C-         ChLIWRK    DIMENSION OF IWRK.  AT LEAST 2*Mcheb + 2.
C-
C-      FAILURE INDICATOR PARAMETER
C-         IFAIL    FAILURE INDICATOR.
C-                     0 - SUCCESSFUL TERMINATION.
C-                     1 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C-                         HAS BEEN VIOLATED -
C-                            Mcheb AT LEAST 1,
C-                            Ncheb = Mcheb + nDers(1) + nDers(2) + ... + nDers(Mcheb),
C-                            ChLWRK AT LEAST 7*Ncheb + 5*IMAX + Mcheb + 2,
C-                            ChLIWRK AT LEAST 2*Mcheb + 2.
C-                     2 - FOR SOME I, nDers(I) IS LESS THAN 0.
C-                     3 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C-                         HAS BEEN VIOLATED -
C-                            XMIN STRICTLY LESS THAN XMAX,
C-                            FOR EACH I, Xch(I) MUST LIE IN THE
C-                               INTERVAL XMIN TO XMAX,
C-                            THE Xch-VALUES MUST ALL BE DISTINCT
C-                     4 - NOT ALL PERFORMANCE INDICES LESS THAN
C-                         ONE, BUT ITMAX ITERATIONS PERFORMED,
C-                     5 - COMPUTATION TERMINATED BECAUSE
C-                         ITERATIONS DIVERGING.
C-
C-
C-      CHECK AND SET ITERATION LIMITS
C-
C-      .. Parameters ..
C-      CHARACTER*6       SRNAME
C-      PARAMETER         (SRNAME='E01AEF')
C-      .. Scalar Arguments ..
C-      DOUBLE PRECISION  XMAX, XMIN
C-      INTEGER           IFAIL, ITMAX, ITMIN, ChLIWRK, ChLWRK, Mcheb, Ncheb
C-      .. Array Arguments ..
C-      DOUBLE PRECISION  chebA(Ncheb), WRK(ChLWRK), Xch(Mcheb), Ych(Ncheb)
C-      INTEGER           nDers(Mcheb), IWRK(ChLIWRK)
      write(*,*)' entering OutputResults_controls_equator_
     *	interpolate_Chebyshev'
C: define for each I from 1 to Mcheb the values of nDers(I), Xch(I),
C                and (Ych(J),J=N+1,N+nDers(I)+1).
C            nDersMAX = MAX(nDersMAX,nDers(I))
C            N = N + nDers(I) + 1
C   *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate_Chebyshev_parChebyshev:'
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
      phia=phiRcylZ(icyl,0)
      enddo
      Mcheb=NrDAT+1
C- too many points
      Mcheb=8
C- try to do in pieces (may be overlapping), a loop over  pieces will be needed
      skipCheb=20
      Mcheb=(NrDAT+1)/skipCheb
      write(*,*)'  NrDAT+1,skipCheb,Mcheb=',NrDAT+1,skipCheb,Mcheb
C-          pause;
      Ncheb=0
      do I=1,Mcheb
      nDers(I)=0
C- non-zero if derivatives are given
      Xch(I)=Rcyld(skipCheb*I-1)
      do J=Ncheb+1,Ncheb+nDers(I)+1
      Ych(J)=phiRcylZ(skipCheb*J-1,0)
      enddo
      Ncheb=Ncheb+nDers(I)+1
      enddo
      XMIN=Xch(1)
      XMAX=Xch(Mcheb)
      IF(Ncheb.LE.ChNMAX.AND.nDersMAX.LE.nDerMX)THEN
      IFAIL=1
      nDersMAX=MAXVAL(nDers)
C- here must be zero
      write(*,*)' nDersMAX=',nDersMAX
C-             call mypause;
      ENDIF
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate_Chebyshev_parChebyshev:'
      write(*,*)' before   IF (Ncheb.LE.ChNMAX  ...)'
      write(*,*)' Ncheb.LE.ChNMAX .AND. nDersMAX.LE.nDersMX',Ncheb,
     *	ChNMAX,nDersMAX,nDersMX
C-          pause;
      IF(Ncheb.LE.ChNMAX.AND.nDersMAX.LE.nDersMX)THEN
      IFAIL=1
      ITMIN=0
      ITMAX=0
C- The next call to E01AEF can be used for obtaining  interpolant q(x) as a series of Chebyshev
C- polynomials with coefficients chebA(i).
C- The polynomial interpolant can subsequently be evaluated for any value of x in the
C- given range by using
C- E02AKF. Chebyshev-series representations of the derivative(s) and integral(s) of q(x) may be
C- obtained by (repeated) use of E02AHF and E02AJF.
      CALL E01AEF(Mcheb,XMIN,XMAX,Xch,Ych,nDers,Ncheb,ITMIN,ITMAX,chebA,
     *	ChWRK,ChLWRK,IWRK,ChLIWRK,IFAIL)
      write(*,*)' after call E01AEF(Mcheb, ...)'
C-             pause;
      WRITE(NOUT,*)
      IF(IFAIL.EQ.0.OR.IFAIL.GE.4)THEN
      WRITE(NOUT,89999)'Total number of interpolating conditions =',
     *	Ncheb
      WRITE(NOUT,*)
      WRITE(NOUT,*)'Interpolating polynomial'
      WRITE(NOUT,*)
      WRITE(NOUT,*)'   I    Chebyshev Coefficient chebA(I+1)'
      DO I=1,Ncheb
      WRITE(NOUT,89998)I-1,chebA(I)
      ENDDO
      WRITE(NOUT,*)
      WRITE(NOUT,*)'  Xch    R   Rth derivative    Residual'
      IY=0
      IRES=nDersMAX+1
      DO I=1,Mcheb
      nDers1=nDers(I)+1
      DO J=1,nDers1
      IY=IY+1
      IRES=IRES+1
      IF(J-1.NE.0)THEN
      WRITE(NOUT,89997)J-1,Ych(IY),ChWRK(IRES)
      ELSE
      WRITE(NOUT,89996)Xch(I),'   0',Ych(IY),ChWRK(IRES)
      ENDIF
      ENDDO
      ENDDO
      ELSE
      WRITE(NOUT,89995)'E01AEF exits with IFAIL =',IFAIL
      write(*,*)' Ifail after Chebyshev E01AEF=',Ifail
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
      write(*,'(a,1p,(5g12.3))')' Xch(I)=',(Xch(I),I=1,MCheb)
C-                call mypause;
C-                pause; -- better than mypause, because exits smoothly on non-enter
      ENDIF
      WRITE(NOUT,89995)'Chebyshev E01AEF exits with IFAIL =',IFAIL
      write(*,*)' Ifail after Chebyshev E01AEF=',Ifail
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
      write(*,'(a,1p,(5g12.3))')' Xch(I)=',(Xch(I),I=1,MCheb)
      do i=1,MCheb
      IFAIL=0
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
C- this computes a function froom Chebyshev expansion:
      CALL E02AKF(NCheb+1,XMIN,XMAX,ChebA,1,NCheb+1,Xch(I),resCheb(I),
     *	IFAIL)
      WRITE(*,'(a,i5,1p,2g12.3)')' ifail Xch(i) testCheb =',ifail,
     *	Xch(i),resCheb(i)
C- corresponds to ptype 6 3 in interactive sm:
      pt(i)=63.d0
      enddo
C-                call mypause;
C-                pause; -- better than mypause, because exits smoothly on non-enter
      ENDIF
C-       CALL E01BAF(N,Xch,Ych,LAMDA,C,LCK,WRK,LWRK,IFAIL);
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate_Chebyshev:'
C: for similar tasks E01BAF, or try the ready 2D spline *
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate:'
      Zcyl=Zd(0)
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
      phia=phiRcylZ(icyl,0)
C-                rho=EOS(phiu*(phiRcylZ(icyl,jcyl)),1)/rho0;
      psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)))
      rho=EOS(phiu*(C-phia-psi),1)
C- in phys. units from EOS
C-            write(*,'(a,1p,4g15.5)')' Rcyl, r, phia, psi=',Rcyl, r(icyl), phia, psi;
C-            write(*,'(a,1p,g15.5)')' rho on equ=',rho;
C-            write(*,'(a,1p,g15.5)')' rho saved =',rhoOn_r_thetaGrid(icyl,JMax);
C-       call mypause;
      enddo
C------- '<--Leaving  Node %_OutputResults_controls_equator:'
C: output CMcyl, CJcyl etc for rho centered *
C------- '-->Entering Node %_OutputResults_controls_rhocentered:'
      CMcyl(0)=0.d0
C- mass within Rcyl(icyl)
      CJcyl(0)=0.d0
C- angular momentum within Rcyl(icyl)
      write(7,'(3a)')'#     Rcyla,              CMcyl(icyl),           
     *	sigma','             sigmaGM,              ERROR','       IFA
     *IL,    CJcyl(icyl)'
      do icyl=1,NrDAT
      Rcyl=Rcyld(icyl)
C- here dp
      if(icyl.GT.1)then
C- now this "if  " is not needed
      Rcylm=Rcyld(icyl-1)
      else
      Rcylm=0.d0
      endif
      Rcyla=half*(Rcyl+Rcylm)
C-            psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)));
      psi=Fpsi(Rcyla,rotlaw(1:len_trim(rotlaw)))
      c1=-Fdpsi_drcyl(Rcyla,rotlaw(1:len_trim(rotlaw)))*Rcyla
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,c1))
      CMcyl(icyl)=CMcyl(icyl-1)
      CJcyl(icyl)=CJcyl(icyl-1)
      sigma=0.d0
      do jcyl=1,NzDAT
      Zcyl=Zd(jcyl)
C- here dp
C
C               if (jcyl==1) then;
C                  write(*,*)' Zcyl=',Zcyl;
C                  call mypause;
C               endif;
C*
      phia=(phiRcylZ(icyl-1,jcyl-1)+phiRcylZ(icyl-1,jcyl)+
     *	phiRcylZ(icyl,jcyl-1)+phiRcylZ(icyl,jcyl))/four
C-                rho=EOS(phiu*(phiRcylZ(icyl,jcyl)),1)/rho0;
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      rhoRcylZ(jcyl,icyl)=rho
C- 1st index for z
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-             call mypause;
      dV=2.d0*(Rcyl**2-Rcylm**2)*dZcyl*pi
C- multiplied by 2 to account for 'southern' hemisphere
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
C-  	       rhodV(jcyl)=2.d0*pi*rho*(Rcyl**2-Rcylm**2);
      	CM=CM+rho*dV
      	sigma=sigma+2.d0*rho*dZcyl
C- multiplied by 2 to account for 'southern' hemisphere
      CMcyl(icyl)=CMcyl(icyl)+rho*dV
      CJcyl(icyl)=CJcyl(icyl)+rho*vphi*Rcyla*dV
      	CJ=CJ+rho*vphi*Rcyla*dV
      	CI=CI+rho*Rcyla**2*dV
C-                P=EOS(phiu*(phiRcylZ(icyl,jcyl)),4)/(rho0*phiu);
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
C-  	       W=W+rho*phiRcylZ(icyl,jcyl)/2d0*dV;
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      enddo
C- loop on z
      Call GillMiller(Zd(1),rhoRcylZ(1,icyl),NzDAT,ANS,ERROR,IFAIL)
      sigmaGM=2.d0*ANS
C- multiplied by 2 to account for 'southern' hemisphere
      sigma2pir(icyl)=2.d0*pi*sigmaGM*Rcyla
      write(7,'(1p,5g20.8,i5,g20.8)')Rcyla,CMcyl(icyl),sigma,sigmaGM,
     *	ERROR,IFAIL,CJcyl(icyl)
      enddo
C- loop on Rcyl
      CallGillMiller(Rcyld(1),sigma2pir(1),NrDAT,CMGM,ERROR,IFAIL)
      write(7,'(a,1p,3g20.8,i5)')'# CM CMGM ERROR IFAIL=',CM,CMGM,
     *	ERROR,IFAIL
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' output controls'
      write(2,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
      write(*,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
      if(Lpolytrope)then
      call polytrope(CN)
      write(2,2)'poly',rho0,r(IB)**delta/alphaPolytrope,Re/
     *	alphaPolytrope,CM*(Re/alphaPolytrope)**3/(4.d0*pi),V*(Re/
     *    alphaPolytrope)**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,
     *    CI*Re**5/(alphaPolytrope**5),
C- *rho0),
     *T/abs(W),Pip/Eth+one,VT
      endif
C-        pause;
C------- '<--Leaving  Node %_OutputResults_controls_rhocentered:'
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
C: output CMcylg, CJcylg etc for rho on cylindrical grid  *
C------- '-->Entering Node %_OutputResults_controls_rhoOnCylGrid:'
      CMcylg(0)=0.d0
C- mass within Rcyl(icyl)
      CJcylg(0)=0.d0
C- angular momentum within Rcyl(icyl)
      write(8,'(3a)')'#     Rcyla,              CMcylg(icyl),    
     *       sigma','             sigmaGM,              ERROR','  
     *     IF AIL,    CJcylg(icyl)'
      Rcylm=0.d0
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
C- here dp
      psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)))
      c1=-Fdpsi_drcyl(Rcyl,rotlaw(1:len_trim(rotlaw)))*Rcyla
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,c1))
      if(icyl.GT.0)then
      CMcylg(icyl)=CMcylg(icyl-1)
      CJcylg(icyl)=CJcylg(icyl-1)
      Rcylm=Rcyld(icyl-1)
      endif
      sigma=0.d0
      do jcyl=0,NzDAT
      Zcyl=Zd(jcyl)
C- here dp
C
C               if (jcyl==1) then;
C                  write(*,*)' Zcyl=',Zcyl;
C                  call mypause;
C               endif;
C*
      phia=phiRcylZ(icyl,jcyl)
C- here phia on cylindrical grid
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      rhoRcylZ(jcyl,icyl)=rho
C- 1st index for z
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-             call mypause;
      dRcyl=Rcyl-Rcylm
      if(jcyl.LT.NzDAT)then
      dV=2.d0*((Rcyl**2-Rcylm**2)+dRcyl**2)*dZcyl*pi
C- multiplied by 2 to account for 'southern' hemisphere
      sigma=sigma+2.d0*rho*dZcyl
C- multiplied by 2 to account for 'southern' hemisphere
      else
      dV=((Rcyl**2-Rcylm**2)+dRcyl**2)*dZcyl*pi
C- not multiplied by 2 since dZcyl is on border
      sigma=sigma+rho*dZcyl
C- not multiplied by 2
      endif
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
C-  	       rhodV(jcyl)=2.d0*pi*rho*(Rcyl**2-Rcylm**2);
      	CM=CM+rho*dV
      CMcylg(icyl)=CMcylg(icyl)+rho*dV
      CJcylg(icyl)=CJcylg(icyl)+rho*vphi*Rcyla*dV
      	CJ=CJ+rho*vphi*Rcyla*dV
      	CI=CI+rho*Rcyla**2*dV
C-                P=EOS(phiu*(phiRcylZ(icyl,jcyl)),4)/(rho0*phiu);
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
C-  	       W=W+rho*phiRcylZ(icyl,jcyl)/2d0*dV;
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      enddo
C- loop on z
      CallGillMiller(Zd,rhoRcylZ(0,icyl),NzDAT+1,ANS,ERROR,IFAIL)
      sigmaGM=2.d0*ANS
C- multiplied by 2 to account for 'southern' hemisphere
      sigma2pir(icyl)=2.d0*pi*sigmaGM*Rcyl
      write(8,'(1p,5g20.8,i5,g20.8)')Rcyl,CMcylg(icyl),sigma,sigmaGM,
     *	ERROR,IFAIL,CJcylg(icyl)
      enddo
C- loop on Rcyl
      Call GillMiller(Rcyld,sigma2pir,NrDAT+1,CMGM,ERROR,IFAIL)
      write(8,'(a,1p,3g20.8,i5)')'# CM CMGM ERROR IFAIL=',CM,CMGM,ERROR,
     *	IFAIL
      do icyl=0,NrDAT
      mcyl(icyl)=CMcylg(icyl)/CMcylg(NrDAT)
      enddo
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' output controls'
      write(8,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *	,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one
     *,VT
      write(*,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *	,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one
     *,VT
      if(Lpolytrope)write(8,2)'poly',rho0,r(IB)**delta/alphaPolytrope,
     *	Re/alphaPolytrope,CM*(Re/alphaPolytrope)**3/(4.d0*pi),V*(Re/al
     *phaPolytrope)**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI*Re**5
     */(alphaPolytrope**5),
C- *rho0),
     *T/abs(W),Pip/Eth+one,VT
C-        pause;
      LinitMcyl=.false.
C- needed for jmscf -- 1st model is rigid, 2nd jmscf
C------- '<--Leaving  Node %_OutputResults_controls_rhoOnCylGrid:'
C------- '<--Leaving  Node %_OutputResults_controls:'
C------- '<--Leaving  Node %_OutputResults:'
      enddo
C------- '<--Leaving  Node %_tableIB:'
      else
C:  -- IB loop decreasing polar radius *
C------- '-->Entering Node %_loopIB:'
      IB=IA
      do while((IB.GT.3))
C  .and. rotlaw=='alpha') .or.
C                    (IB>0.6*IA .and. rotlaw=='rigid'));*
C- assuming IB<=9999
      IF(.NOT.(IB.LT.10))GOTO 09987
      write(str,'(a,I1)')'000',IB
      GOTO 09984
09987 CONTINUE
      IF(.NOT.(IB.LT.100))GOTO 09986
      write(str,'(a,I2)')'00',IB
      GOTO 09984
09986 CONTINUE
      IF(.NOT.(IB.LT.1000))GOTO 09985
      write(str,'(a,I3)')'0',IB
      GOTO 09984
09985 CONTINUE
      write(str,'(I4)')IB
09984 CONTINUE
      open(26+IB,file='phi'//str(1:len_trim(str))//'.res')
C:  iterate phi etc. for lower IB, which means
C                             a flatter configuration, so faster rotation *
C------- '-->Entering Node %_FindPhi:'
      write(*,*)' entering _FindPhi'
C-            call mypause;
      LogInit=.true.
      write(*,'(a,1p,g15.5,i7)')'rho0, IB',rho0,IB
      write(*,'(a,1p,3g15.5)')'phi(IA,JMax),phi(I0,JMax) =',phi(IA,JMax)
     *	,phi(I0,JMax)
C-           'phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)))=',phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)));
C- write of Fpsi is forbidden if there is write in it
      write(*,'(2a)')' rotlaw=',rotlaw(1:len_trim(rotlaw))
      C=(phi(IA,JMax)+Fpsi(one,rotlaw(1:len_trim(rotlaw)))-phi(I0,JMax)
     *	*EOS(zero,3)/EOS(rho0,3)
C- EOS(*,3) means enthalpy H
     *)/(one-EOS(zero,3)/EOS(rho0,3))
C- eq.17 in AA290(1994)674
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      write(*,'(a,1p,3g15.5)')'phiu,EOS(rho0,3),C=',phiu,EOS(rho0,3),C
C-            call mypause;
      n_max=300
      epsMain=1.d-12
C- to data file
      LogIter=.true.
      nStep=1
      do while(LogIter.and.(nStep.LE.n_max))
      C0=C
      Callnew_phi(
C-r,alpha,theta,dtheta, -- now in module cylinderGrid
     *C,LogInit,LogRapid,rotlaw(1:len_trim(rotlaw)))
C- last argument is the rotation law
C              if(C>C0)then;
C                C=C0+min(abs(C-C0),0.03d0*abs(C0));
C              else;
C                C=C0-min(abs(C-C0),0.03d0*abs(C0));
C              endif;*
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      if(abs((C0-C)/C).LT.epsMain)LogIter=.false.
      write(*,*)' nStep,C,phi(I0,JMax):',nStep,C,phi(I0,JMax)
      nStep=nStep+1
      enddo
C------- '<--Leaving  Node %_FindPhi:'
      if(IB.le.0.2*IA.or.rotlaw.EQ.'rigid')then
      IB=IB-nint(0.01*IA)
      else
      if(IB.le.0.4*IA)then
      IB=IB-2*nint(0.01*IA)
      else
      if(IB.le.0.6*IA)then
      IB=IB-4*nint(0.01*IA)
      else
      IB=IB-8*nint(0.01*IA)
      endif
      endif
      endif
C        do i0=0,IMax0;
C                 rcyl=float(i0)/IMax0;
C                 im=powrsmart(rcyl,one/delta)/dr;
C                 ip=im+1;
C                 rcylm=im*dr;
C                 rcylp=ip*dr;
C                 xm=(rcylp-powrsmart(rcyl,one/delta))/dr;
C                 xp=(powrsmart(rcyl,one/delta)-rcylm)/dr;
C                 phiOut=xm*phi(im,JMax)+xp*phi(ip,JMax);
C                 write(26+IB,'(1p,3e15.4)')xm, xp, phiOut;
C               enddo;
C             *
      enddo
C------- '<--Leaving  Node %_loopIB:'
      endif
C-**      open(16,file='phi');
C-**      write(16,*) phi0;
C-**      close(16);
      stop
C:  *
1     Format(i5,4(1pe10.2),1pe28.20)
2     Format(a,1p,2g17.8,e17.8,g17.8,12g15.6,I5,g15.6)
22    Format(1p,2g17.8,e17.8,g17.8,13g15.6,I5,g15.6)
C-
99999 FORMAT(1X,8F9.4)
99998 FORMAT(F5.2)
99997 FORMAT(1X,I16,F12.4,I11,F12.4)
99996 FORMAT(1X,I39,F12.4)
89999 FORMAT(1X,A,I4)
89998 FORMAT(1X,I4,1p,g20.3)
89997 FORMAT(5X,I4,1p,g12.1,g17.6)
89996 FORMAT(1X,F4.1,A,1p,g12.1,g17.6)
89995 FORMAT(1X,A,I2,A)
      end
C:  SUBROUTINE GillMiller(X,Y,N,ANS,ER,IFAIL)
C              for numerical integration of table data *
      SUBROUTINE GillMiller(X,Y,N,ANS,ERROR,IFAIL)
      DOUBLEPRECISION ANS,ERROR
      INTEGER IFAIL,N
C-     .. Array Arguments ..
      DOUBLEPRECISION X(N),Y(N)
      IFAIL=1
      CALL D01GAF(X,Y,N,ANS,ERROR,IFAIL)
      GOTO(09983,09982,09981,09980),IFAIL+1
      stop' D01GAF produced wrong IFAIL'
      GOTO 09979
09983 CONTINUE
      if(abs(ANS).GT.0.d0)then
      if(abs(ERROR/ANS).GT.0.1d0)WRITE(*,'(1X,A,1p,G12.4,A,G12.4)')'
     *	Integral = ',ANS,'     Estimated error = ',ERROR
      endif
      GOTO 09979
09982 CONTINUE
      WRITE(*,*)'Less than 4 points supplied'
      GOTO 09979
09981 CONTINUE
      WRITE(*,*)'Points not in increasing or decreasing order'
      GOTO 09979
09980 CONTINUE
      WRITE(*,*)'Points not all distinct'
09979 CONTINUE
      return
      end
C:  *
      SUBROUTINE D01GAF(X,Y,N,ANS,ER,IFAIL)
C-
C-     THIS SUBROUTINE INTEGRATES A FUNCTION (Y) SPECIFIED
C-     NUMERICALLY AT N POINTS (X), WHERE N IS AT LEAST 4,
C-     OVER THE RANGE X(1) TO X(N).  THE POINTS NEED NOT BE
C-     EQUALLY SPACED, BUT SHOULD BE DISTINCT AND IN ASCENDING
C-     OR DESCENDING ORDER.  AN ERROR ESTIMATE IS RETURNED.
C-     THE METHOD IS DUE TO GILL AND MILLER.
C-
C-     NAG COPYRIGHT 1975
C-     .. Scalar Arguments ..
      DOUBLEPRECISION ANS,ER
      INTEGER IFAIL,N
C-     .. Array Arguments ..
      DOUBLEPRECISION X(N),Y(N)
C-     .. Local Scalars ..
      DOUBLEPRECISION C,D1,D2,D3,H1,H2,H3,H4,R1,R2,R3,R4,S
      INTEGER I,NN
C-     .. Executable Statements ..
      ANS=0.0D0
      ER=0.0D0
      IF(N.GE.4)GOTO20
      IFAIL=1
      RETURN
C-
C-     CHECK POINTS ARE STRICTLY INCREASING OR DECREASING
C-
20    H2=X(2)-X(1)
      DO80I=3,N
      H3=X(I)-X(I-1)
      IF(H2*H3)40,60,80
40    IFAIL=2
      RETURN
60    IFAIL=3
      RETURN
80    CONTINUE
C-
C-     INTEGRATE OVER INITIAL INTERVAL
C-
      D3=(Y(2)-Y(1))/H2
      H3=X(3)-X(2)
      D1=(Y(3)-Y(2))/H3
      H1=H2+H3
      D2=(D1-D3)/H1
      H4=X(4)-X(3)
      R1=(Y(4)-Y(3))/H4
      R2=(R1-D1)/(H4+H3)
      H1=H1+H4
      R3=(R2-D2)/H1
      ANS=H2*(Y(1)+H2*(D3/2.0D0-H2*(D2/6.0D0-(H2+2.0D0*H3)*R3/12.0D0)))
      S=-(H2**3)*(H2*(3.0D0*H2+5.0D0*H4)+10.0D0*H3*H1)/60.0D0
      R4=0.0D0
C-       write(*,'(a,1pe12.3)')' init ans:',ans;
C-
C-     INTEGRATE OVER CENTRAL PORTION OF RANGE
C-
      NN=N-1
      DO120I=3,NN
      ANS=ANS+H3*((Y(I)+Y(I-1))/2.0D0-H3*H3*(D2+R2+(H2-H4)*R3)/12.0D0)
C-          write(*,'(a,i5,1pe12.3)')' i ans:',i,ans;
      C=H3**3*(2.0D0*H3*H3+5.0D0*(H3*(H4+H2)+2.0D0*H4*H2))/120.0D0
      ER=ER+(C+S)*R4
      IF(I.NE.3)S=C
      IF(I.EQ.3)S=S+2.0D0*C
      IF(I-N+1)100,140,100
100   H1=H2
      H2=H3
      H3=H4
      D1=R1
      D2=R2
      D3=R3
      H4=X(I+2)-X(I+1)
      R1=(Y(I+2)-Y(I+1))/H4
      R4=H4+H3
      R2=(R1-D1)/R4
      R4=R4+H2
      R3=(R2-D2)/R4
      R4=R4+H1
      R4=(R3-D3)/R4
120   CONTINUE
C-
C-     INTEGRATE OVER FINAL INTERVAL
C-
140   CONTINUE
      ANS=ANS+H4*(Y(N)-H4*(R1/2.0D0+H4*(R2/6.0D0+(2.0D0*H3+H4)*
     *	R3/12.0D0)))
C-       write(*,'(a,1pe12.3)')' final ans:',ans;
      ER=ER-H4**3*R4*(H4*(3.0D0*H4+5.0D0*H2)+10.0D0*H3*(H2+H3+H4))/
     *	60.0D0+S*R4
      ANS=ANS+ER
      IFAIL=0
      RETURN
      END
